P5176 公约数

首先有一个结论:

gcd(ij,ik,jk)=gcd(i,j)gcd(j,k)gcd(i,k)gcd(i,j,k)gcd(i\cdot j,i\cdot k,j\cdot k)=\frac{gcd(i , j )\cdot gcd(j , k) \cdot gcd(i,k)}{gcd(i,j,k)}

证明如下:

先将 i,j,ki,j,k 用唯一分解定理展开的次方分别为 w1,w2,w3w_1,w_2,w_3

gcd(ij,jk,ik)=i=1kpimin(w1i+w2i,w2i+w3i,w1i+w3i)gcd(i \cdot j,j \cdot k,i \cdot k)=\sum_{i=1}^k {p_i}^{min(w_{1i}+w_{2i},w_{2i}+w_{3i},w_{1i}+w_{3i})}

gcd(i,j)gcd(j,k)gcd(i,k)=i=1kpimin(w1i,w2i)+min(w2i,w3i)+min(w1i,w3i)gcd(i,j) \cdot gcd(j,k) \cdot gcd(i,k)=\sum_{i=1}^k {p_i}^{min(w_{1i},w_{2i})+min(w_{2i},w_{3i})+min(w_{1i},w_{3i})}

gcd(i,j,k)=i=1kpimin(w1i,w2i,w3i)gcd(i,j,k)=\sum_{i=1}^k {p_i}^{min(w_{1i},w_{2i},w_{3i})}


gcd(ij,ik,jk)=gcd(i,j)gcd(j,k)gcd(i,k)gcd(i,j,k)gcd(i\cdot j,i\cdot k,j\cdot k)=\frac{gcd(i , j )\cdot gcd(j , k) \cdot gcd(i,k)}{gcd(i,j,k)}

i=1kpimin(w1i+w2i,w2i+w3i,w1i+w3i)=i=1kpimin(w1i,w2i)+min(w2i,w3i)+min(w1i,w3i)i=1kpimin(w1i,w2i,w3i)\Leftarrow \sum_{i=1}^k {p_i}^{min(w_{1i}+w_{2i},w_{2i}+w_{3i},w_{1i}+w_{3i})}=\frac{\sum_{i=1}^k {p_i}^{min(w_{1i},w_{2i})+min(w_{2i},w_{3i})+min(w_{1i},w_{3i})}}{\sum_{i=1}^k {p_i}^{min(w_{1i},w_{2i},w_{3i})}}

i=1kmin(w1i+w2i,w2i+w3i,w1i+w3i)=i=1kmin(w1i,w2i)+min(w2i,w3i)+min(w1i,w3i)min(w1i,w2i,w3i)\Leftarrow \sum_{i=1}^k {min(w_{1i}+w_{2i},w_{2i}+w_{3i},w_{1i}+w_{3i})}=\sum_{i=1}^k {min(w_{1i},w_{2i})+min(w_{2i},w_{3i})+min(w_{1i},w_{3i})-min(w_{1i},w_{2i},w_{3i})}

这是一个轮换式 , 不妨令 w1iw2iw3iw_{1i} \le w_{2i} \le w_{3i} , 则有

i=1kw1i+w2i=i=1kw1i+w2i+w1iw1i\Leftarrow \sum_{i=1}^k {w_{1i}+w_{2i}}=\sum_{i=1}^k {w_{1i}+w_{2i}+w_{1i}-w_{1i}}

i=1kw1i+w2i=i=1kw1i+w2i\Leftarrow \sum_{i=1}^k {w_{1i}+w_{2i}}=\sum_{i=1}^k {w_{1i}+w_{2i}}

所以原结论成立。


i=1nj=1mk=1pgcd(ij,ik,jk)×gcd(i,j,k)×(gcd(i,j)gcd(i,k)×gcd(j,k)+gcd(i,k)gcd(i,j)×gcd(j,k)+gcd(j,k)gcd(i,j)×gcd(i,k))\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p\gcd(i\cdot j,i\cdot k,j\cdot k)\times \gcd(i,j,k)\times \left(\frac{\gcd(i,j)}{\gcd(i,k)\times \gcd(j,k)}+\frac{\gcd(i,k)}{\gcd(i,j)\times \gcd(j,k)}+\frac{\gcd(j,k)}{\gcd(i,j)\times \gcd(i,k)}\right)

i=1nj=1mk=1pgcd(i,j)gcd(j,k)gcd(i,k)gcd(i,j,k)×gcd(i,j,k)×(gcd(i,j)gcd(i,k)×gcd(j,k)+gcd(i,k)gcd(i,j)×gcd(j,k)+gcd(j,k)gcd(i,j)×gcd(i,k))\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p\frac{gcd(i , j )\cdot gcd(j , k) \cdot gcd(i,k)}{gcd(i,j,k)}\times \gcd(i,j,k)\times \left(\frac{\gcd(i,j)}{\gcd(i,k)\times \gcd(j,k)}+\frac{\gcd(i,k)}{\gcd(i,j)\times \gcd(j,k)}+\frac{\gcd(j,k)}{\gcd(i,j)\times \gcd(i,k)}\right)

i=1nj=1mk=1pgcd(i,j)gcd(j,k)gcd(i,k)×(gcd(i,j)gcd(i,k)×gcd(j,k)+gcd(i,k)gcd(i,j)×gcd(j,k)+gcd(j,k)gcd(i,j)×gcd(i,k))\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p gcd(i , j)\cdot gcd(j , k) \cdot gcd(i,k)\times \left(\frac{\gcd(i,j)}{\gcd(i,k)\times \gcd(j,k)}+\frac{\gcd(i,k)}{\gcd(i,j)\times \gcd(j,k)}+\frac{\gcd(j,k)}{\gcd(i,j)\times \gcd(i,k)}\right)

i=1nj=1mk=1pgcd(i,j)2gcd(j,k)2gcd(i,k)2\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p gcd(i , j)^2 \cdot gcd(j , k)^2 \cdot gcd(i,k)^2

p×i=1nj=1mgcd(i,j)2+n×i=1mj=1pgcd(i,j)2+m×i=1nj=1pgcd(i,j)2p \times \sum_{i=1}^n\sum_{j=1}^m gcd(i,j)^2+n \times \sum_{i=1}^m\sum_{j=1}^p gcd(i,j)^2+m \times \sum_{i=1}^n\sum_{j=1}^p gcd(i,j)^2

现在只需要处理出 i=1nj=1mgcd(i,j)2\sum_{i=1}^n\sum_{j=1}^m gcd(i,j)^2 ,就可以用三次反演解决。

s=min(n,m)s = min(n,m)

i=1nj=1mgcd(i,j)2\sum_{i=1}^n\sum_{j=1}^m gcd(i,j)^2

k=1sk2i=1nj=1m[gcd(i,j)=k]\sum_{k=1}^{s}k^2\sum_{i=1}^n\sum_{j=1}^m [gcd(i,j)=k]

k=1sk2i=1nkj=1mk[gcd(i,j)=1]\sum_{k=1}^{s}k^2\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{k} \rfloor} [gcd(i,j)=1]

k=1sk2d=1skμ(d)nkdmkd\sum_{k=1}^{s}k^2 \sum_{d=1}^{\lfloor \frac{s}{k} \rfloor } \mu(d) \lfloor \frac{n}{kd} \rfloor \lfloor \frac{m}{kd} \rfloor

k=1sd=1skk2μ(d)nkdmkd\sum_{k=1}^{s}\sum_{d=1}^{\lfloor \frac{s}{k} \rfloor } k ^2 \mu(d) \lfloor \frac{n}{kd} \rfloor \lfloor \frac{m}{kd} \rfloor

k=1skd=1sk2μ(d)nkdmkd\sum_{k=1}^{s}\sum_{kd=1}^{s} k ^2 \mu(d) \lfloor \frac{n}{kd} \rfloor \lfloor \frac{m}{kd} \rfloor

T=kdT=kd ,

T=1skTk2μ(Tk)nTmT\sum_{T=1}^{s}\sum_{k|T} k^2 \mu(\frac{T}{k}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor

T=1sf(T)nTmT\sum_{T=1}^{s}f(T)\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor

其中 f(n)=knk2μ(nk)f(n)=\sum_{k|n} k^2 \mu(\frac{n}{k}) , 由卷积的性质得该函数为积性函数,考虑如何用线性筛求出其前缀和。

如果 pp 是一个质数,那么有:

f(p)=12×μ(p)+p2×μ(1)=p21f(p)=1^2 \times \mu(p)+p^2 \times \mu(1)=p^2-1

f(pk)=i=0kp2i×μ(pki)f(p^k)=\sum_{i=0}^k p^{2i} \times \mu(p^{k-i})

又因为当 ki2k-i \ge 2 时,μ(pki)=0\mu(p^{k-i})=0,所以所有有贡献的 i>k2i > k - 2 , 即 i=ki=ki=k1i=k-1

f(pk)=p2kp2(k1)=(p21)p2(k1)=f(p)p2(k1)f(p^k)=p^{2k}-p^{2(k-1)}=(p^2-1)p^{2(k-1)}=f(p)p^{2(k-1)}

观察发现括号内 pp 减少了 k1k-1 次方,系数增加 p2(k1)p^{2(k-1)}

所以 f(pk)=f(pk1)×p2f(p^k)=f(p^{k-1}) \times p^2

整理一下得:

f(i×prime[j])={f(i)×f(prime[j])(i,prime[j])=1f(i)×prime[j]2(i,prime[j])=1f(i \times prime[j])=\begin{cases} f(i) \times f(prime[j]) & (i,prime[j])=1 \\\\ f(i) \times prime[j]^2 & (i,prime[j])\not=1 \end{cases}

复杂度Θ(n+tn)\Theta(n+t\sqrt n)

#include <cstdio>
#include <iostream>
using namespace std; 
#define Mod 1000000007

template<typename _T>
void Read( _T &x ) {
	x = 0; int f = 1;
	char s = getchar( );
	for( ; s < '0' || s > '9' ; s = getchar( ) ) f = s == '-' ? -f : f;
	for( ; s >= '0' && s <= '9' ; s = getchar( ) ) x = x * 10 + s - '0';
	x *= f;
}
template<typename _T>
void Write( _T x ) {
	if( x < 0 ) putchar('-') , x = -x;
	if( x >= 10 ) Write( x / 10 );
	putchar( x % 10 + '0' );
}

const int MAXN = 20000000;
int t , n , m , p , k , prime[ MAXN + 5 ] , f[ MAXN + 5 ];
bool vis[ MAXN + 5 ];

void sieve( ) {
    f[ 1 ] = 1;
    for( int i = 2 ; i <= MAXN ; i ++ ) {
        if( !vis[ i ] ) {
            prime[ ++ k ] = i;
            f[ i ] = ( 1ll * i * i - 1 ) % Mod;
        }
        for( int j = 1 ; j <= k && 1ll * i * prime[ j ] <= MAXN ; j ++ ) {
            vis[ i * prime[ j ] ] = 1;
            if( i % prime[ j ] == 0 ) {
                f[ i * prime[ j ] ] = 1ll * f[ i ] * prime[ j ] % Mod * prime[ j ] % Mod;
                break;
            }
            f[ i * prime[ j ] ] = 1ll * f[ i ] * f[ prime[ j ] ] % Mod;
        }
    }
    for( int i = 2 ; i <= MAXN ; i ++ )
        f[ i ] = ( f[ i ] + f[ i - 1 ] ) % Mod;
}

int solve( int n , int m ) {
    int d = min( n , m ) , Ans = 0;
    for( int l = 1 , r ; l <= d ; l = r + 1 ) {
        r = min( n / ( n / l ) , m / ( m / l ) );
        Ans = ( Ans + 1ll * ( f[ r ] - f[ l - 1 ] + Mod ) % Mod * ( n / l ) % Mod * ( m / l ) % Mod ) % Mod;
    }
    return Ans;
}

int main( ) {
    sieve( );
    Read( t );
    while( t -- ) {
        Read( n ) , Read( m ) , Read( p );
        int Ans = ( 1ll * n * solve( m , p ) % Mod + 1ll * m * solve( n , p ) % Mod + 1ll * p * solve( n , m ) % Mod ) % Mod;
        Write( Ans ) , putchar('\n');
    }
    return 0;
}